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The long-term stability of extrasolar system HD 37124. Numerical 
study of resonance effects. 
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ABSTRACT 

We describe numerical tools for the stability analysis of extrasolar planetary systems. In par- 
ticular, we consider the relative Poincare variables and symplectic integration of the equations 
of motion. We apply the tangent map to derive a numerically efficient algorithm of the fast 
indicator MEGNO (a measure of the maximal Lyapunov exponent) that helps to distinguish 
chaotic and regular configurations. The results concerning the three-planet extrasolar system 
HD 37124 are presented and discussed. The best fit solutions found in earlier works are stud- 
ied more closely. The system involves Jovian planets with similar masses. The orbits have 
moderate eccentricities, nevertheless the best fit solutions are found in dynamically active 
region of the phase space. The long term stability of the system is determined by a net of 
low-order two-body and three-body mean motion resonances. In particular, the three-body 
resonances may induce strong chaos that leads to self-destruction of the system after Myrs 
of apparently stable and bounded evolution. In such a case, numerically efficient dynamical 
maps are useful to resolve the fine structure of the phase space and to identify the sources of 
unstable behavior. 

Key words: extrasolar planets — ^Doppler technique — stars individual HD 37124 — N-body 
problem — numerical methods 



1 INTRODUCTION 

Understanding the extrasolar planetary systems has became a ma- 
jor challenge for contemporary astronomy. One of the most dif- 
ficult problems in this field concerns the orbital stability of such 
systems. Usually, the investigations of long-term evolution are the 
domain of direct, numerical integrations. The stability of extraso- 
lar systems is often understood in terms of the Lagrange definition 
implying that orbits remain well bounded over an arbitrarily long 
time. Other definiti ons may be f ormulated as well, like the astro- 
nomical stability ( Lissauerll 19991) requiring that t he system persists 
over a very long, Gyr time-scale, or Hill stability dSzebehelvll 19841) 
that requires the constant ordering of the planets. In our studies, we 
prefer a more formal and stringent approach relat ed to the funda- 
mental Kolmogorov-Arnold-Theorem (KAM), see lArnold JT973) . 
Planetary systems, involving a dominant mass of the parent star and 
significantly smaller planetary masses, are well modeled by close- 
to-integrable, Hamiltonian dynamical systems. It is well known, 
that their evolution may be quasi-periodic (with a discrete number 
of fundamental frequencies, forever stable), periodic (or resonant; 
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stable or unstable) or chaotic (with a continuous spectrum of fre- 
quencies, and unstable). In the last case, initially close phase trajec- 
tories diverge exponentially, i.e., their Maximum Lyapunov Char- 
acteristic Exponent (MLCE, denoted also with cr) is positive. In 
general, the distinction between regular and chaotic trajectories is a 
very difficult task that may be resolved only with numerical meth- 
ods relying on efficient and accurate integrators of the equations of 
motion. 

The main motivation of this paper is to describe numerical 
tools that are useful for studies of th e dynamical stab ility and to 
apply them to the HD 37124 system jVogt et alJboOSi) . We recall 
the fundamentals of relative canonical Poincare variables as - in 
our opinion - one of the best frameworks for symplectic integra- 
tors, jnies^^anonic^v^iables are well suited for the construction 
of a iLaskar & Robutell j200lh composition m ethod that improves 
a clas sical Wisdom-Holman (W-H) algorithm l lWisdom & HolmM] 
Il99lh . We supplement the integrator with a propagator of the as- 
sociated symplectic tan gent ma p that approximates the solution of 
variational equations tMikkola & Innanenll 19991 ). Finally, we com- 
pare two fast indicators that reveal the character of phase trajecto- 
ries. The first one is a relatively simple method for resolving funda- 
mental frequencies and spectral properties of a close-to-integrable 
H amiltonian system - a so called Spect ral Number (SN), invented 
bv lMichtchenko & Ferraz- The second indicator be- 

longs to the realm of the Lyapunov expo nent based algorithms; w e 
chose the numerical tool developed by ICincotta & Simd ( l2000h ; 
ICincotta et afl ^20031) under the name of MEGNO. In this work. 
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we refine the algorithm of MEGNO that ma kes explicit use of the 
symplectic tangent map l lGozdziewskill2003l) . 

As a non-trivial application of the presented numerical tools, 
we consider the 3-planet system hosted by the HD 37124 star 
dVogt et all |2003) . It has been discovered by the radial veloc- 
ity (RV) technique. The recent model of the RV observations of 
HD 37124 predicts three equal Jovian type planets with masses 
~ 0.6 mj in orbits with moderate eccentricities. In such a case, 
the application of symplectic integrators without regularization is 
particularly advantageous thanks to the numerical efficiency (long 
time-steps) and accuracy (the total energy does not have a secular 
error and the angular momentum integral is conserved). The num- 
ber of multi-planet systems resembling the architecture of the So- 
lar system increaseo Hence, our approach may be useful in other 
cases. 



'H(P,P): 
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where the mutual distance A,j is 



Mlp;-P,i 



(2) 
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and k stands for the Gaussian gravity constant. 

The Hamiltonian ^ admits six integrals of barycenter 



N 

J]m,p, = 0, Y/i = ^, (4) 

!=0 (=0 

as well as the energy integral "K = const, and three angular momen- 
tum integrals 



:£p,XP,-: 



const. 



(5) 



2 NUMERICAL TOOLS 

According to the classical results of celestial mechanics, the A'- 
body problem has only 10 integrals of motion for all N > 2; they 
consist of 6 integrals of barycenter, 3 integrals of angular momen- 
tum and the energy integral. The integrals of barycenter play a very 
particular role in the studies of an A'-body system dynamics. First, 
they define the origin of an inertial reference frame in terms of the 
mutual distances and velocities of the bodies considered, thus dis- 
missing the need of some extrinsic absolute frame. But what is 
more important, being linear forms of coordinates and momenta 
they allow a unique reduction of the system, lowering the number 
of degrees of freedom by three, with no loss of information. This is 
why we can solve the relative two-body problem and then recover 
the motion of both masses with respect to their center of mass. And 
this is why we can approximately solve the heliocentric motion of 
planets, recovering the barycentric evolution a posteriori. 

Within the framework of Hamiltonian mechanics, the reduc- 
tion is usually achieved by means of a transformation to one of the 
two common variable typ es: relative Jacob i variables, or "heliocen- 
tric" Poincare variables jWhittakeill 19521 Ch. XIII). We focus on 
the latter set, because it offers the best choice in many aspects. We 
introduce the basic ideas related to the Poincare variables and we 
derive them as a Mathieu transformation; this way is simpler and 
more intuitive than the procedure b ased u pon a gener ating function 
that w as presented by IWhittakeil ( Il952l p. 343) or iDuncan et al.l 
( Il998l) . Then we discuss the setup of the Hamiltonian within the 
framework of Wisdom-Holman type integrators. 

2.1 Poincare variables basics 

Let us consider a system consisting of Af -I- 1 material points with 
masses mo, . . . , m^r. We define a barycentric position vector p 6 R-"^ 
and its canonical conjugate momentum P 6 R-''^ as 



(1) 



Then, the barycentric equations of m otion can be derived from the 
Hamiltonian function ( lLaskai|[l990h 
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The integrals are usually exploited as the accuracy control tool 
when the differential equations of motion 



5P ' 



P = - 



(6) 



are solved numerically. 

Instead of solving the 6A' + 6 order system ((SJ, it is often de- 
sirable to study the relative motion of A' bodies with respect to 
the material point mo. But if the integration method to be applied 
is symplectic, it is necessary to use the Hamiltonian equations of 
motion, hence the necessity of defining a canonical transformation 
(p,P,'K) ^ (r,R,'7C) with new, relative coordinates r and momenta 
R. 

A naive, straightforward approach would consist in postulat- 
ing r,- = p,- - po for all This leads to ro = and the Jacobian matrix 
of the transformation becomes singular; such a transformation can- 
not be canonical regardless of the choice of the momenta. However, 
the di fficulty can be ea sily circumvented if we change the definition 
of rn . IPoincarel i 1 89d) proposed 



fori= 1, 



(7) 



r/ = Pz-Po 

ro = Po, 

whereas iDuncan et alj ( Il998h chose ro as the position of the 
barycenter in some arbitrary inertial frame. For the barycentric 
system, the latter choice amounts to a "differentiable zero" ro = 

(Z/Io ™') 2jIo "^'P' ~ turns out, that both starting points lead 
to the same result, so we continue our presentation assuming the 
Poincare choice Q. 

Retaining ro as the position of the reference body mo with 
respect to the barycenter, we can perform a canonical extension 
of the time-independent point tran sformation Q, requesting the 
Mathieu transformation condition ( IWhittakerlll95l p. 301) 



P dp = R dr, 

or, explicitly. 



(8) 



N 



^ Pi • dp/ = Ro -dpo + Yj • ('IP'- - 'IPo) • (9) 

;=o !=i 

Equating the coefficient of each differential dp/ to zero, we find the 
new momenta 

R/ = P/, for ! = !,..., A', 



Ro = 2f=oP/ = 



(10) 



Planets Encyclopedia, |http://exoplanets.eu| 



The fact, that Ro = is a direct consequence of the integrals of 
barycenter (|4j. 
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Equations Q and l llOt . that are due to lPoincar^ ( Il896h . define 
the canonical relative variables. The coordinates r consist of the 
positions with respect to the reference body niQ, save for the ro 
that is measured with respect to the barycenter. The momenta R 
are measured with respect to the barycenter, save for Rq that can 
be understood as measured with respect to otq and hence it is zero 
(although with nonvanishing partials with respect to P,). 

The most important feature is that the new Hamiltonian "K, 
obtained by the simple substitution of the transformation equations 
into "H, does not depend neither o n the coordinates, nor on the mo- 
ment a of mo. Indeed, one obtains ('Poinca re|[T893 : IWhittakeJl 1 9521 : 
Hagihara 1970; Deprit 1983; Laskar 1990 019911) 



1^/1 1\ 2 1 ^ ^ 
2-f— '\mo mil niQ ^ Ar' 
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-z 



k niQini 



N N 
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wher e A,-/ = ||p/ -p,i| = 
ll998l ; IChamb"erslll999l) 

I ^ I ^ 
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or, equivalently jPuncan et al.l 
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An important fact to be remembered is that the Hamiltonian 'K 
has the form i ll lb or ( 112b only if the substitution of the barycenter 
integrals Q has been performed. Thus it cannot serve to obtain the 
equations for Vq or Ry. However, once we know all the remaining 
r; and R,-, the Tq values can be easily computed from Equations {4j 
and 0, whereas Rq = by the definition. For all the remaining 
bodies 



(13) 



and we will assume that i = l,...,N throughout the rest of this pa- 
per. 

A remarkable property of the Poincare variables is 



G = J]p,xP, = J]r,xR;, 



(14) 



1=0 



1=1 



which means, that the total angular momentum of the reduced sys- 
tem of N bodies evaluated by means of the Poincare variables is 
the same as the angular momentum of A' + 1 bodies evaluated in the 
barycentric frame. 

If the reference body mass itiq » m,-, the Hamiltonian can 
be easily partitioned into the unperturbed, Keplerian part and a 
small perturbation proportional to the greatest of m, . From the point 
of view of analy tical theories using these variables jYuasa & Horil 
ll979l ; lHorill985l) , it is preferable to split 'K into the sum 



7C = ^^'"+^J'", 

where jLaska]ill990lll99lh 

2^ mam: ' ^ 
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(15) 



(16) 
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The principal part TCq"' defines A' relative two-body problems. The 

perturbation TCj"' involves not only the mutual interactions be- 
tween the minor bodies m,-, but also the momenta related terms 
that replace the usual "indirect part" of the perturbing function 
prese nt in noncanonical relative (A'+ l)-body problem jPoincard 
Il905l) . It is due to this term, that the Poincare variables were con- 
sidered somehow handicapped; the objection that velocities are no 
longer tangent to the momenta became almost a proverb, although 
many non-inertial reference frames have the same property, the 
restricted three-body problem being the best example. This ob- 
jection has fortunately ceased to be taken seriously; for example, 
|Ferraz-Mello et al. ( 2004) successfully use orbital elements evalu- 
ated from the Poincare momenta R. In this paper we will use alter- 
natively two types of orbital elements. The osculating elements are 
computed by the usual two body formulae from astrocentric posi- 
tions r/ and velocities r/; we use them in all plot labels and orbital 
data tables. But in the calculation of spectral numbers or in the def- 
initions of resonance arguments, we use orbital elements computed 
from r,- an d mom, /(mo -Hm/)R/, calling them contact elements after 
iBrumber^ ( Il99 ll ). [In fact, the transformation between astrocen- 
tric positions - barcentric momenta and the co ntact elements ca n 
be still done with the usual two-body formulae ( iMorbidellil2002l) 1. 
The former are commonly used in literature, whereas the latter olfer 
a better behavior from the dynamical point of view. The superiority 
of contact elements results from the fact, that the reference frame 
of Cartesian momenta is inertial, hence the influence of noninertial 
forces is reduced to purely kinematical contribution. The inferred 
Keplerian angles (w, SI, M) and the conjugat e momenta can be in- 
terpreted as canonical Delaunay's elements ( lMorbidelIill2003) , so 
the derivation and interpreation of the fundamental frequencies is 
straightforward. 



2.2 Symplectic integration: two is a company 

With the advent of symplectic integrators ba sed on the Wisdom- 
Holman approach jWisdom & Holmanlll99lh . the Poincare vari- 
ables became an at tractive framework for the numerical studies of 
planetary systems touncanetal]|l998l ; IChambers|[T99^ . First of 
all, similarly to the lacobian coordinates, they reduce the number 
of equations of motion by 6 with respect to the barycentric prob- 
lem. Thanks to the possibility of splitting '7C into the main part and 
a perturbation, they also allow the construction of a symplectic in- 
tegrator with the local truncation error proportional to the product 
of nii/mo, hence a larger integration step h can be applied. 

Given a Hamiltonian M = Mo +eMi with a small parame- 
ter E, the W-H integrator is based on the alternating application of 
maps <l)() ,-(r,R) and Oi ,-(r, R) that represent the solutions of equa- 
tions of motion derived from Mo and eMi alone, on the interval 
from to to to + T. Moreover, each of the Hamiltonian parts should 
admit (a possibly simple) analytical solution of the equations of 
motion. 

In the numerical applications, Hamilto nian TC is typical! 
split differently than in Equation ( |15b , namely jPuncan et al.ll99l 
IChamberslll999h 



n< = 'Ko+n<u 



where 



'Ko 



iylR2 y 

/=i (=1 



k momj 



(18) 



(19) 
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'Kx = — 
2mo 



I N \ 
Vi=l 



N N ,2 

/=i j=i+i 



(20) 



The unperturbed part TCo has now a different meaning: it still leads 
to A' relative two-body problems 

' flR/ in, 



R, = - 



dr, 



(21) 
(22) 



but this time they are the restricted two-body problems with neg- 
ligible masses m; or a fixed center of gravity: 



k^rriQ 



(23) 



The perturbation part "Ki remains proportional to mi/mo 
and it is still a func tion of b oth coordinates r and momenta R. 
iDuncan et al. I ( Il998l) and then IChambersI ( 1 19991) considered it an 
obstacle, so they further split TCi into 



TCi ='7Cii(R) + '7Ci2(r), 

obtaining elementary "kick" maps <I>ii_t 



(24) 



where r, stands for r,(/o), r' stands for r,(/o + r), and similarly for 

r/r; 

Separating '/C according to Equations l |18l l, l |19l l, and i20\ of- 
fers a possibility of using only two maps for a W-H integrator. Al- 
though our partitioning seems to be more compact and clear than 
the "heliocentric-democratic" scheme, we note that both mappings 
are practically equivalent as far as the CPU cost is concerned. Yet 
the rule "two is a company, three is a crowd" holds true in the realm 
of symplec tic integrators for perturbed systems: according to the 
theorem of ISuzukj ( fl99lh any symplectic composition method of 
an order higher than 2 must necessarily involve stages with negative 
sub-steps that amplify accumulation of roundoff errors. This holds 
true for a composition of maps derived f rom split ting the H amilto- 
nian into any number of terms. However. iLaskar & Robutell l l200lh 
found a family of methods designed for two-terms perturbed sys- 
tems with TC = ^ + eS where the negative sub-steps are avoided. 
Their integrators do not contradict the results of Suzuki: formally 
they remain second order methods, but in contrast to other W-H 
methods with local truncation errors O(et^), their errors have a 
form 0(e^t^ + et"). So, for sufficiently small perturbation e, the 
Laskar-Robutel methods may behave like higher order integrators 
in certain domain of stepsize t although no negative substeps were 
introduced. 



r 



mo pi 



r; = R, 

and Oi2,T- 



r; 



R, 



N 



k mimj 



(r,-r,). 



(25) 
(26) 

(27) 
(28) 



In the formulas of both maps we add a prime to the symbols stand- 
ing for the values of coordinates and momenta at to + t, whereas 
unprimed symbols refer to the values at tg. 

As the effect of the partition I l24t . the classical "leapfrog" 

*T~*l,T/2°*0,T°*l,T/2, (29) 

was replaced by 

1>r ~ l,T/2 ° f 12,T/2 ° *0,T ° ® 12,T/2 ° ® 1 1,t/2- (30) 

According to lDuncan et alj i 19981), the ordering of On and (S>i2^t 
is insignificant, and, indeed, [Chamber j fT999l) interchanged them, 
using 

fr ~<I>12,T/2°f ll,T/2°fo,T°®ll,T/2°<I'l2,T/2- (31) 

The interchange of the maps is justified by the fact that "Ki \ and 
'Ki2 commute, i.e. the Poisson bracket CTCn ; 'K12] = 0. In these cir- 
cumstances 



*1.T = ®12,T°®11.T = ®11.T°®12,T, 



(32) 



and we can concatenate both maps obtaining 9i j in a compact 
form. 



mo 



N 



R' 



R/ - k^ mi T E ■ 



(33) 



(34) 



2.3 Tangent maps 

Whenever a differential correction of initial conditions or the com- 
putation of sensitivity indicators is required, the use of tangent 
maps becomes indispensable. Keplerian map Oq and its associate 
ta ngent map can be computed according to a comprehensive recipe 
bv iMikkola & InnanenI ( [l999l) . The propagation of a tangent vector 



dr 
<5R 



(35) 



under the action of the "kick" map Oi amounts to multiplying it by 
the Jacobian matrix DO]. Resulting expressions are simple: 



Sr' 

5r; 

where 



(5r+ — y 5R,-, 



(r;-r,)-<5,- 



6ij + - 



3<5A/y 



5ij = 6rj-6ri. 



(36) 



(37) 



(38) 



2.4 Angular momentum integral 

It can be easily demonstrated, that any composition of maps Oq and 
O] conserves the angular momentum integral l ll4t . This property is 
guaranteed by the conservation of G by each map separately. Re- 
calling that Oo and ®l define exact solutions of motion generated 
by "Ko and 'TCi respectively, we can simply check that 

{G;7Co| = {G;7Ci| = 0. (39) 

The proof of Eq. l |39l l is straightforward. Starting from the defini- 
tion of G, we use the linearity of Poisson brackets and the Leibnitz 
identity to write 

N 

iG;7Col = E iR';'^ol-RiXir,;'7Co| = 
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^ 5r,- dRi 



(40) 



Then we substitute the right-hand sides of Equations ( 121b and l l22b . 
concluding that all vector products vanish and indeed |G; 'TCq} = 0. 
A similar procedure demonstrates jG; TCi } = 0. 

Of course, from practical point of view the conservation of G 
is only up to computer roundolf errors. 



2.5 Chaoticity indicators 

To detect unstable motions in the phase space, many numerical 
tools are available. Concerning the dynamics of close-to-integrable 
Hamiltonian systems, they can be roughly divided in two classes: 
spectral algorithms t hat resolve the fundamental frequencies and/or 
their diffusion rates ( Laskar 1993: Sidlichovs kv & Nesvomvll 19971 : 
iMichtchenko & Ferraz-Mello 2001) . and methods based on the di- 
vergence rate of initially close p hase trajectories, expressed in 
terms of the Lyapunov exponents jSenettin et al.|[l980l : iFroeschld 
Il984h . 

In this work, among the the spectral tools, we ch oose the 
method invented by IMichtchenko & Ferraz- Melld ( l200lh : its idea 
is genuinely simple — to detect chaotic behavior one counts the 
number of frequencies in the FFT-spectrum of an appropriately 
chosen dynamical signal. We deal with conservative Hamiltonian 
systems; so in a regular case, the spectrum of fundamental frequen- 
cies is discrete and we obtain only a few dominant peaks in the 
FFT spectrum. Chaotic signals do not have well defined frequen- 
cies, and their FFT spectrum is very complex. The number of peaks 
in the spectrum above some noise level p (typically, p is set to a few 
percent of the dominant amplitude) tell us on the character (regular, 
chaotic) of the of the s ystem. 

The method by IMichtchenko & Ferraz-Mellol ( 1200 ih does 
not have as strong theor etical f oundations as the Frequency 
Map Analysis (FMA) b y iLaskai] ( | 1993,) or the Fourier Modi- 
fied Transform (FMT) by lSidlichovskv & Nesvoni^ ( Il997h which 
are considered as rigorous and efficient tools. We did some 
comparative tests of the later algorithm with MEGNO already 
('Oozdziewski & Konacki 2004). Here, we choose the method of 
[Michtchenko & Ferraz-Mello C2001,) for its appealing simplicity 
and because we used it in the former papers devoted to the analysis 
of the RV data. In that way, we can compare the results directly. 
In our code, the spectral signals analyzed with the FFT are related 
to canonical Poincare elements, so the fundamental frequencies are 
well defined. Moreover, we resolve the chaotic and regular signals 
by comparing the number of significant peaks in the FFT spectrum, 
thus a very precise determination of the fundamental frequencies is 
not critical. Actually, in this work we also show that the algorithm 
has same drawbacks and should be applied with care. 

The basic tool to discover exponentially unstable bounded or- 
bits, i.e. chaotic orbits, is the Maximum Lyapunov Characteristic 
Exponent (MLCE) cr. Numerical symplectic integration methods 
are fixed step algorithms, so we can restrict our discussion to the 
iterations of a discrete map = '^'"(q, that generates a sequence 
of state vectors ^„ consisting of coordinates and their conjugate 
momenta. The direct computation of the MLCE is based on the 
analysis of the tangent vectors 6„ that evolve under the action of a 
linear tangent map dn = (D<t>)"do. Asymptotically, the MLCE value 
is given by 



limiyinf^ 

\Sk-i 



If cr converges to some positive value, we conclude that the nominal 
orbit and some initially close orbit diverge exponentially at the 
rate exp((Tf). Two practical difficulties arise when the direct defini- 
tion ( I41t is used: the convergence of cr is often very slow, and it is 
difficult to tell how small should be the final value of a to consider 
it cr = 0. 

A large variety of methods has been proposed to overcome 
the problem of slowly convergent MLCE estimates. The authors 
prefer the so called MEGNO (Mean E xponential Growth factor of 
Nearby Orbits) indicator proposed by ICincotta & Sim3 ( |2000|) - 
that choice is justified by th e successful application of this method 
in ou r previous works (e.g.. lBreiter et al.ll2005l : lGozdziewski et al.l 
'2006, and the re ferences therein). Th e definition of MEGNO for a 
discrete map is ( ICincotta et al.ll2003h 



Y{n)=-yy{k), 
n i—i 



k=\ 



where 



,(«) = -VHn(-^ 



(42) 



(43) 



If the iterates of the discrete map refer to the moments of time 
separated by the stepsize h, the discrete map MEGNO function Y{n) 
asymptotically tends to 

Y„ = ahn + b, 

with a = 0,b = 2 for a quasi-periodic orbit, a = b = foi a stable, 
isochronous periodic orbit, and a = |cr, b = for a chaotic orbit. 
Thus we can indirectly estimate the MLCE on a finite time interval, 
but the weight function k in the sum ( |43! ) reduces the contribution of 
the initial part of the tangent vector evolution, when the exponential 
divergence is to small to be ob served behind other linear and non- 
linear effects ( lMorbidellill200 2l). Thus, fitting the straight line to the 
final part of F(w), we obtain good estimates of cr from a relatively 
shorter piece of trajectory than in the direct MLCE evaluation. 

In practical application, one can use a more con venient form 
of Eqs. (|42j and ^ proposed bv lBreiter et"al] j2005l) 



Y(n) 



(n-l)Y(n-l)+y(n) 



y(n) = —y(n- l) + 2ln(^ 



(44) 
(45) 



(41) 



with the initial setup ^(O) = Y(0) = 0. The fact that only the ratio 
6nlSn-\ is significant, as well as the linearity of tangent map, al- 
lows to avoid the overflow of 6n thanks to occasional normalization 
of the tangent vector length to (5„ = 1 performed after the ratio of 
dn/S„-\ has been evaluated. 



3 STABILITY OF THE HD 37124 PLANETARY SYSTEM 

As a non-trivial application of the presented algorithms and the 
illustration of difficulties arising in the dynamical analysis of the 
long-term stability of multiplanet configurations, we choose the 
HD 37124 extrasolar syste m. The discovery o f two Jovian plan- 
ets has bee n announced by iButler et"al] ( 1200 ih and confirmed by 
IVogt et al] (2005). At first, t he system seemed to be well modeled 
by a 2-planet configuration dButler et alj|200ll) . However, new ob- 
servations lead to two-planet fits with ~ 0.7 and a catastrophi- 
cally unst able configuration . Moreover, with the updated RV ob- 
servations, IVogt etaljjioo^) found much better model of 3 planets 
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Table 1. The bet-fit astro-centric, osculating Keplerian elements of a 
stable HD 37124 planetary configuration at the epoch of the first obser- 
vation (o=JD2,451,0420.047. Mass of the parent star is 0.78 mp The 
fit has been refined with GAMP over ~ 5 • 10-*Pd- See iVoet et al'llioo^; 
iGozdziewski et alj2006h and Fig.|2]for more details. 



Parameter 


planet b 


planet c 


planet d 


msini [mj] 


0.62447 


0.56760 


0.71194 


a [AU] 


0.51866 


1.61117 


3.14451 


e 


0.07932 


0.15267 


0.29775 


w [deg] 


138.405 


268.863 


269.494 


M{to) [deg] 


259.011 


109.545 


124.113 






0.938 




Vo[m s-'] 




7.629 




rms [m s"'] 




3.39 





with similar masses of ~ 0.6 mj in low-eccentric orbits. The best 
fits have the rms ~ 4 m/s, in agreement with the internal accuracy 
of the data. However, the best-fit orbital solution, both in the the 
kinematic Keplerian model, and in more realistic W-body simula- 
tion (see Table 1), lies close to the collision line of planets c and 
d. Note, that we define the collision line in terms of semi-axes and 
eccentricities as adl + e^) = aj(l - ej). This line marks the zone in 
which the mutual interactions of massive companions can quickly 
destabilize the system. 

How to interpret the RV measurements remains an open ques- 
tion. The dynamical long-term stability of the planetary system is 
the most natural requirement of a configuration consistent with ob- 
servations. Yet the three-planet model is parameterized by at least 
16 parameters, even assuming that the system is coplanar. For that 
reason the search for the best fits fulfilling the constraints of sta- 
bility is a difficult task. It can be resolved in different ways. For 
example, we may try to find dynamically stable solutions in the 
vicinity of the formal best fit configurations (the latter are often 
unstable). However, examining the stability of configurations in 
that neighborhood, we have no reasons to expect that the stable 
fits are optimal in the statistical sense. Another approach relies on 
the elimination of unstable fits during the fitting process, through 
penalizing unstable solutions with a large val ue of (Xy)^^"^. This 
method, described in IGozdziewski et alj j2006h . is dubbed GAMP 
(Genetic Algorithm with MEGNO Penalty). It was shown, that 
such an approach is particularly useful in modeling resonant or 
close-to-resonant planetary configurations. Unfortunately, the al- 
gorithm cannot give definite answer when we want to resolve the 
ixiy^^ shape of strictly stable solutions in detail. The penalty term 
in the (XyY^^ function relies on a signature of the system stabil- 
ity, expressed through the fast indicator. Due to significant CPU 
overhead, the fast indicator in the minimizing code can be only cal- 
culated over relatively short time, typically lO-' orbital periods of 
the outermost planet. Moreover, the code can converge to unstable 
best fits that appear stable on that short time scale. Hence, at the 
end of the search, we have to examine the stability of the individual 
best fits in the obtained ensemble of solutions, over the time-scale 
of relevant mean motion and secular resonances. 



3.1 Long-term stability of the best-fit configurations 

Th e GAMP a nalysi s of the RV dat a of HD 37124 publishe d 
bv lVogt etal.l j2005l) was presented in IGozdziewski etaU j2006l) . 
About 100 of best fit solutions were found yielding ixlV^^ < I.I, 
the rms ~ 4 m/s, and stable in the sense that their MEGNO signa- 



tures are close to 2 up to ~ 1000 - 2000 orbital periods of the out- 
ermost planet. Due to heavy CPU requirements, the time-span to 
resolve MEGNO in the GAMP code cannot be set very long. The 
short integration time ~ 10^ allows only to eliminate strongly 
chaotic configurations, typically leading to collisions between plan- 
ets and/or with the parent star. The best fit solutions were found 
in a dynamically active region of the phase space, spanned by a 
number of low-order mean motion resonances (MMRs) between 
the two outermost Jovian companions, like 5c:2d, 8c:3d, or llc:4d 
(see Fig.|2|and dynamical maps presented in Figs. l3l4l5l6t . In par- 
ticular, close to the collision lines, the low-order MMRs overlap, 
giving rise to the global instability zone. 

Yet we should be aware that two-body MMRs with char- 
acteristic time scale ~ 10^-10^ orbital periods of the outermost 
planet are not the only source of instability in the multi-planet sys- 
tem. Already when we deal with three-planet configurations, the 
strong instabilities may be ge nerated by three-body MMRs or by 
long-term secular r e sonance s iNesvom^ & Morbidellilll99g,ll999l : 
iMurrav et alj|l998l : IGuzzoI|2006|) . In such instance, appropriately 
longer integration time is necessary to detect the unstable solutions. 

This issue is illustrated in Fig. [T] We choose one of the best 
fits with initial osculating, astrocentric Keplerian elements at the 
epoch of the first observation in terms of tuples {m [mj], a [AU], 
e, (jj [deg], M [deg]): (0.593, 0.519, 0.0058, 303.360, 95.060), 
(0.558, 1.615, O.IOI, 315.621, 70.279), and (0.690, 3.193, 0.26111, 
255.848, 142.886) for planets b, c, and d, respectively. These ini- 
tial conditions had (Xv)^^^ ~ 0.98 and an rms about 4 m/s. In the 
time range covered by the GAMP integrations (and up to ~ 10^ Pj, 
i.e. ~ 60,000 yr), the configuration appears strictly regular because 
the indicator quickly converges to 2 (the top-left panel in Fig[T|l. 
Nevertheless, after the transient time, the MEGNO starts to grow 
linearly at the rate of cr/2 ~ 2- lO^^yr"' (where cr is the MLCE 
of the solution, see bottom-left panel in Fig. [T](. Actually, after a 
relatively long time ~ 15 Myr, the chaotic motion leads to a colli- 
sion between planets c and d (the top-right panel) due to a sudden 
increase of both eccentricities up to 0.6. The elimination of such so- 
lutions during an extensive GAMP-like search on a Myrs interval 
would be very difficult. 

Looking for the source of such dramatically unstable behavior, 
we perform the frequency analys is of the orbits with the MET by 
ISidlichovskv & Nesvomvl jl997h . Denoting the proper mean mo- 
tions by nb>"c, and nd> respectively, we found that 

;jb-8nc +Vnd a; -0.4° /yr, 

clearly indicating the three-body MMR of the first order, and we 
label it with -I- lb : -8c : +7d. The time evolution of the critical ar- 
gument 6 = A\,- 8/ic + 7/ici is illustrated in the bottom-right panel of 
Fig. El The circulation of the critical angle alternates with libration, 
indicating the separatrix crossings that explain chaotic evolution. 

The presented example has inspired us to follow a two-stage 
procedure in modeling the RV data. First, with a GAMP-like code 
we look for many best-fit solutions, ideally, approximating the 
global shape of (Xy)^^^ and simultaneously stable, at least over a 
relatively short period of time. At that stage the stability constraints 
cannot be tight, not only due to significant CPU requirements but 
also because we should not discard weakly chaotic solutions. Such 
configurations may be bounded over very long time, longer by or- 
ders of magnitude than their Lyapunov time Tl = I /cr. In the next 
step, we either refine the search in a zone bounded by the previously 
found fits with much longer integration times (still numerically ex- 
pensive), or we examine each fit with long-term direct integrations 
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Figure 1. Evolution of the HD 37124 system selected best fit related to the theree-body MMR (initial conditions listed in the text). Left: evolution of MEGNO 
over a short (top) and a long (bottom) time. The straight line is the least-square fit to Y{t) = (cr/2)t + b. Top right: contact eccentricities during ~ 15 Myr. 
Bottom right: the critical argument of the three-body MMR +lb : -8c : -^7d. 



and/or evaluate a fast indicator signature, like the MLCE, Spectral 
Number, or the diffusion of fundamental frequencies. 

Here, for each solution with ixlY^'^ < 1-1> we computed its 
MEGNO signature. The integration time span is about of 37 Myr 
- long enough to detect the relevant chaotic three-body reso- 
nances and strong secular resonances. Here, and in the experi- 
ments des cribed later on, we use the SBAB3 integrator scheme by 
iLaskar & R obutel (2001). The time-step is 4 days. The secular pe- 
riods in the given range of Cd are quite short, ~ 10** yr, nevertheless 
we can expect that dynamical effects of potentially active secular 
resonances could be detected after thousands of such characteristic 
periods, hence counted in 10^-10^ yr. Figure |2] illustrates the re- 
sults. The quality of fits in terms of ixiY^'^ is marked by the size 
of circles (better fits have larger circles). Red (medium grey) cir- 
cles are for stable, quasi-periodic solutions. In that case the system 
may be stable over a very long time. Blue (dark grey) circles are 
for chaotic solutions that led to collisions between planets and that 
did not survive during the integration time (the integrations are in- 
terrupted if any of the eccentricities increases above 0.66). Finally, 
small yellow (light grey) circles mark all configurations (not neces- 
sarily regular) that survived, remaining bounded during the maxi- 
mal integration time. Clearly, most of solutions with initial > 0.2 
are both chaotic and unbounded. Nevertheless, some chaotic solu- 
tions appear on the borders of stable regions as well. Generally, the 
distribution of fits gives us a clear image of the border of global 
instability of the system, relatively far from the collision zone. 



3.2 Fine structure of the phase space 

Figure [3] compares the sensitivity of MEGNO and the Spectral 
Number when we use the same integration time, ~ 10^ yr ^ 1.6- 
\(f Pi. In the case of the SN map, we did the EFT on A' = 2'*^ steps 
of 64 d, counting the number of spectral lines above 1% of the 
largest amplitude in the signal of /(f) = a^{t)esu^iA^{t), where 
and /Ic denote the contact semi-major axis and mean longitude of 
planet c. 



Both dynamical maps present the same region of the phase 
space, in the neighborhood of the best fit. Note, that this particular 
solution has been refined with GAMP integration over time 5 ~ 
10"* P(] that is about of 2 orders of magnitude longer than in the 
set of selected solutions. The resolution of the maps is the same: 
480 X 120 data points; the map coordinates are usual astrocentric 
osculating Keplerian elements. Most of best-fits from Fig.[2]lie in 
the region covered by Fig. [3] 

Both maps reveal a number of unstable resonances. Yet the 
SN map involves some artifacts (mo/re-like patterns) related to a 
low value of the noise level parameter p. Within the same inte- 
gration time, the MEGNO map reveals relatively more fine details 
than SN. In particular, the sophisticated border of the collision zone 
appears to be more sharp and shifted towards smaller ej. We can 
also find some fine resonance lines entirely absent in the SN map. 
For instance, there is a fine structure on the right-hand side of the 
8c : 3d MMR. In order to investigate that instability, we choose 
the initial condition marked with small crossed circle and labeled 
with a. The results of the MET frequency analysis of this solution 
tell us that the structure is related to the +2b : - 12c : -l-3c MMR. The 
time evolution of the related critical argument is illustrated Fig. [8^. 

Next, we computed close-ups of the dynamical map within the 
rectangle labeled I in Fig. [3] These maps are shown in Fig.|4] This 
time we increased pXo 5% and the total number of steps has been 
doubled {N = 2^") in order to avoid the moire artifacts. But once 
again the "concurrent" MEGNO map (the right panel of Fig. ^ 
calculated over the same total time seems to offer a better represen- 
tation of the phase space. Interestingly, the best fit data (Table I) 
seem to lie on the border of a chaotic zone spanned by many over- 
lapping resonances. A close-up of that area, marked with rectan- 
gle II in Fig.|4l is shown in Fig.|5] Clearly, even a very small change 
of parameters of the outermost planet may push the system into a 
strongly chaotic state. It also illustrates the good performance of 
the GAMP algorithm that was able to locate and preserve the fit in 
an extremely narrow island of stable motions. 

A closer inspection of the area III in Fig. |4] reveals a multi- 
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2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 



Figure 2. The long-term stability of the best-fits obtained in the GAMP search bv lGozdziewski et al.l )2006l) . The osculating elements at the epoch of the first 
observation (JD 2,450,420.047) are projected onto the (flrf,erf)-plane (i.e., the semi-major axis vs the eccenfiicity of the outermost planet). The quality of fits in 
terms of (xl)^^^ < 1.1 and rms< 4 m/s is marked by the size of blue (dark grey) and red (medium grey) circles (larger circles have better fit quality). The best 
fit, self-consistent Newtonian configuration obtained without stability constraints, in terms of quintuples (m [mj], a [AU], e, at [deg], M [deg]) at the epoch of 
the first observation, is (0.619, 0.519, 0.088, 141.91, 257.34), (0.565, 1.663, 0.104, 331.88, 67.71), and (0.732, 2.947, 0.378, 283.33, 95.32) for planets b, c, 
and d, respectively, with velocity offset 7.53 m/s. The blue circles are for chaotic solutions that did not survive the integration time of ~ 37 Myr. The red circles 
are for regular' solutions — in that case the MEGNO converged to 2. The small yellow (light grey) circles mark configurations that survived the integration. 
The best stable fit found is marked with an arrow, its osculating elements are given in Table 1. Some dominant MMRs of planets c and d are marked with 
dashed vertical lines, according to the third Kepler law, and labeled accordingly. The red curve marks the collision line of the two outermost orbits. 



O.S 1 1.5 2 2.5 3 1.8 2 2.2 2.4 2.6 2.8 3 




Figure 3. Dynamical maps of the Spectral Number (left) and the MEGNO indicator (right) computed in the neighborhood of the best-fit solution to the RV 
data of HD 37124 with the integration time ~ 10^ yr. The elements of the best fit (see Table 1) are labeled with the crossed circle. The rectangle (I) marks the 
borders of the close-up shown in Fig.|4] The stabiHty of orbits is color-coded: in both maps, yellow (pale grey) means strongly chaotic and unstable solutions; 
regular configurations are marked black in the SN-map and dark blue (dark grey) in the MEGNO map. 



tude of weakly unstable solutions. To show such structures in more 
detail, we computed a close-up of that area (Fig. |6j. The resolu- 
tion of the maps is 200 x 200 data points, the integration interval is 
~ 3 ■ 10^ yr a; 5 • 10^ Pd- The time step at the left column is 16 days 
and it provides the relative error of the total energy at the level of 
10"^ . Apparently, in spite of relatively short integration time, the 
map uncovers sophisticated structure of the two-body and three- 
body resonances. To identify some of them, we choose initial con- 
dition marked by small crossed circles and labeled in the map with 
b, c, and d, respectively. 

For initial condition b we found that «(, - 11 /ic + I5n^ ~ 
-0.1°/yr, i.e., indicating three-body MMR + lb : -11c : + I5d; for 
initial condition c we have got \0n^-21n^ » -0.3° /yr, correspond- 
ing to the -I- 10c : -27d MMR of the outer giants; and for initial 
condition d we have nb -"c - 12nd ~ -0.2° /yr, indicating the three- 



body + lb:-lc:-12d MMR, respectively. All these resonances ex- 
cite chaotic configurations. That can be demonstrated through ob- 
serving the time evolution of their critical angles (see Fig[8j5,c,d). 
In all those instances, we found that the libration alternates with 
circulation of these angles, hence confirming that the relevant con- 
figurations are close to the resonance separatrices. 

A particularly interesting star-like structure can be be seen 
around the initial condition d (the upper-left panel in Fig. |6]l. In 
that area, the two-body 10c : 27d MMR and many week three-body 
resonances are active, for instance, -l-lb;-llc:-l-15dw -1.35°/yr, 
-t-lb ; -Ic : -12d « -0.2°/yr, +2h : -12c : H-3d « -1.5°/yr, -F3b : 
-13c : -9d a; -1.7°/yr, and -I- 10c : -27d ~ 1.2°/yr. 

One might be tempted t o attribut e this sophisticated structure 
to the so called Arnold web ( ICincotta |[2002). Indeed, a closer look 
at the branches of the web shows new fine details and extremely 



© 0000 RAS, MNRAS 000, 000-000 



The long-term stability of extrasolar system HD 37124. Numerical study of resonance effects. 9 



O.S 1 l.S 2' S5 3 1.8 2 2.2 2.4 2.6' 'Z.^ S 



11:4 
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logSN 
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3.12 3.13 3.14 3.15 3.16 [AU] 3.12 3.13 3.14 3.15 3.16 [AU] 

Figure 4. Dynamical maps in terms of the Spectral Number (the left panel) and the MEGNO indicator (the right panel) computed in the region marked with 
small rectangle I in Fig. |4] The integration time ~ 3 X 10^ yr is equivalent to ~ 4.8 • 10'' orbital periods of the outermost planet. The resolution of the maps 
500 X 120 data points. The stability of orbits is color-coded, see the caption to the previous figure. 
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Figure 5. Dynamical maps in terms of the MEGNO indicator computed 
in the region marked by small rectangle in Fig. \4\ The integration time is 
~ 3 X lO' yr that is equivalent to ~ 4.8 ■ 10* orbital periods of the putative 
outermost planet. The resolution of the plot is 200 X 320 data points. 



complex dynamical structure, in that zone, illustrated in the close- 
up map around initial condition d, Fig. |7] But the truth is that this 
particular structure is mostly spurious, and it occurred due to an 
improper choice of the integration step. To shed more light on that 
issue, we show the map of the relative error of the total energy (the 
left-bottom panel in Fig.|6]l. The coincidence of higher energy error 
streaks (bottom left) with instability patterns detected by MEGNO 
(top left) is not conclusive by itself, but when we recompute both 
maps using a smaller time step of 10 days (panels in the right col- 
umn of Fig.|6j, we notice that the web patterns of higher MEGNO 
disappear (Fig. (6] top right) and the energy error map significantly 
flattens (Fig.|6] bottom right). We conclude that two additional res- 
onance lines that crossed a t d were generated by the so-called 'step- 
size resonances' jRauch & Hohnan 1999) between proper frequen- 
cies of the system and the sampling frequency of the constant step 
integrator. The effect of step-size resonance in a constant step inte- 
grator can be avoided either by using a sufficiently small integration 
step or by the application of high-order schemes. Unfortunately, 
both approaches lead to more time consuming algorithms. 

As we could expect, the symplectic scheme outperforms the 
classical integration algorithms. For instance, we found the the 
ME GNO code driven by th e Bulirsh-Gragg-Stoer ODEX integra- 
tor (|HaireL&^^ne3[l99l) with the relative accuracy set to 10 
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Figure 7. Close-up's of the MEGNO dynamical map shown in the left 
panel of Fig.|6]The integration time ~ 2x lO' yr is equivalent to ~ 3.2- 10* 
orbital periods of the putative outermost planet d. The resolution is 200 X 
480 points. The MEGNO is computed by the symplectic algorithm with the 
time-step of 16 days. 



requires a similar CPU time, as the Laskar-Robutel SBAB3 scheme 
with 4 days step-size, but the former leads to a much larger, secu- 
larly growing energy error (which is larger by 2-3 orders of magni- 
tude). 



4 CONCLUSIONS 

The use of Poincare variables in the studies of the dynamics of 
close-to integrable planetary systems offers many advantages. The 
variables are canonical and offer a simple form of a reduced Hamil- 
tonian. The Hamiltonian can be split into a sum of two separately 
integrable parts: the Keplerian term and a small perturbation. As 
such, it can serve to construct a symplectic integrator based on any 
m odem composition metho d, including the recent ones invented 
bv'Laskar & Robutell ( I2OO1I) . The tangent map computed with the 
same integration scheme provides an efficient way of computing the 
estimate of maximal Lyapunov exponent in terms of relatively re- 
cent fast indicator MEGNO. The method proves to be much more 
efficient than general purpose integrators (like the Bulirsh-Stoer- 
Gragg method). Besides, it provides the conservation of the inte- 
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Figure 1: 

Figure 6. Close-up's of the MEGNO dynamical map shown in Fig. O within rectangle labeled with III, illustrating the fine structure of the phase space. The 
integration time ~ 3 X lO' yr is equivalent to ~ 4.8 ■ lO"* orbital periods of the putative outermost planet d. The resolution is 200 X 200 points. Panels in the top 
are for the MEGNO computed by the symplectic algorithm: the left panel is for the time-step of 16 days, the right panel if for the time step of 10 days. The 
bottom row is for the relatitve error of the total energy, for the same time steps, respectively. 



a) e = 2^^ -nX^ + -I- 05^, -1-3 C3^ -I- 3 b) e = A-t, - 1 1>.^ -H uX^ -3 m^^-m^ - 




r[10'yr] r[10'yr] 



Figure 8. Time-evolution of the critical arguments of some resonances illustrated in the dynamical maps in Figs.|4]and|6] The panels are labeled accordingly 
with initial conditions marked by small crossed circles in these dynamical maps. See the text for more details. 
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grals of energy and the angular momentum that is crucial for re- 
solving the fine structure of the phase space. From the practical 
point of view, the symplectic algorithms are relatively simple for 
numerical implementation. 

Using the numerical tools, we investigate the long term sta- 
bility of extrasolar planetary system hosted by HD 37124. The or- 
bital p^anietersjnjhesetof our best, self-consistent Newtonian 
fits I Gozdz iewski et al. 2006') are in accord with the discovery pa- 
per dVogt et al. 2005). Nevertheless, the observational window of 
the system is still narrow and the derivation of the model consis- 
tent with observations is difficult and, in fact, uncertain. The dy- 
namical maps reveal that the relevant region of the phase space, 
in the neighborhood of the mathematically best fit, is a strongly 
chaotic and unstable zone. The fitting algorithm (GAMP) that re- 
lies on eliminating strongly unstable fits founds solutions with a 
similar quality [in terms of Of^)'^^] that yields the formal solu- 
tion. Moreover, they are shifted towards larger semi-major axes 
and much smaller eccentricities of the outermost planet. The orbital 
evolution of two outer planets is confined to a zone spanned by a 
number of low-order two-body and three-body MMRs. In particu- 
lar, the three-body MMRs may induce very unstable behaviors that 
manifest themselves after many Myrs of an apparently stable and 
bounded evolution. To deal with such a problem, the stability of the 
best fits should be examined over a time-scale that is much longer 
than the one required when only the two-body MMRs are consid- 
ered. In accord with the dynamical maps, the stable fits to the RV of 
HD 37124 should have small eccentricity of the outermost planet d, 
not larger than 0.2-0.3. Moreover, the stable configurations of the 
HD 37124 system are puzzling. The best-fit mathematical three- 
planet model is surprisingly distant, in the phase space of initial 
conditions, from the zone of stable solutions consistent with the 
RV. It remains possible that other bodies are present in the system 
and the three-planet model is not adequate to explain the RV vari- 
ability, in spite that it provides apparently perfect fits. Yet, isolated 
initial conditions or even sets of best-fit solutions do not provide 
a complete answer on the system configuration. Then the fast in- 
dicator approach is essential and helpful to resolve the dynamical 
structure of the phase space. 

The results of our experiments confirm and warn that all nu- 
merical methods should be applied with great care. All symplec- 
tic methods are constant step integrators. In that case one should 
be cautious about the possibility of generating spurious resonance 
webs. A proper way to avoid them is to repeat computations with a 
different integration step in order to detect step-dependent patterns. 
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